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We present and discuss a general density-matrix description of energy-dissipation and decoherence 
phenomena in open quantum systems, able to overcome the intrinsic limitations of the conventional 
Markov approximation. In particular, the proposed alternative adiabatic scheme does not threaten 
positivity at any time. The key idea of our approach rests in the temporal symmetrization and 
coarse graining of the scattering term in the Liouville-von Neumann equation, before applying the 
reduction procedure over the environment degrees of freedom. The resulting dynamics is genuinely 
Lindblad-like and recovers the Fermi's golden rule features in the semiclassical limit. Applications 
, to the prototypical case of a semiconductor quantum dot exposed to incoherent phonon excitation 

' peaked around a central mode are discussed, highlighting the success of our formalism with respect 

to the critical issues of the conventional Markov limit. 
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INTRODUCTION 



In spite of the quantum-mechanical nature of the electron and photon dynamics in the core region of typical solid- 
state nanodevices, the overall behavior of such quantum systems is often characterized by a non-trivial interplay 
between phase coherence and energy relaxation/ dephasing phenomenai^i^ This fundamental topic is getting more 
and more critical in the last years, due to the continuous downscaling in size that state-of-the-art technology allows 
for. The proper modelling of, e.g., charge-carrier dynamics in cascade nanodevices,-^ rather than of spin/charge- 
based implementations of quantum-computation algorithms in solid-state logic gateSf^ are just a couple of examples 
demanding for a unitary microscopic description. Not to forget the still unsolved issue represented by the effects due 
to the presence of spatial boundaries (contacts),^ 

Contrarily to purely atomic and/or photonic systems -where decoherence phenomena may be successfully described 
via phenomenological adiabatic-decoupling procedures- the quantitative modelling of realistic solid state devices 
requires to account for both coherent and incoherent processes, on equal footing. To this aim, motivated by the 
power and flexibility of the semiclassical kinetic theory^ in describing a large variety of interaction mechanisms, a 
quantum generalization of the standard Boltzmann collision operator has been proposedji the latter, obtained via 
the conventional Markov (CM) limit, describes the evolution of the reduced density matrix in terms of in- and out- 
scattering super-operators. However, contrary to the semiclassical case, such collision super-operator does not preserve 
the positive-definite character of the density-matrix operator because of its non-temporal symmetric construction. 

This serious limitation was originally noted by Spohn and co-workers^ three decades ago; in particular, they clearly 
pointed out that the choice of the adiabatic decouplin g st rategy is definitely not unique, and only one among the 
diverse possibilities, developed by Davies in Ref.s could be shown to preserve positivity: it was the case 

of a "small" subsystem of interest interacting with a thermal environment, and selected through a partial trace 
reduction. Unfortunately, the theory was restricted to finite-dimensional subsystems only (i.e., iV- level atoms), and 
to the particular projection scheme of the partial trace. 

Inspired by the pioneering papers by Davies and co-workers, aim of the present work is to propose an alternative and 
more general adiabatic procedure. In particular, we shall further explore the path recently proposed in Ref.s (T2|.(T3| 
which, in the discrete-spectrum case reduces to Davies' model,^- for diagonal states gives the well known Fermi's 
golden rulCj^^^'^'^ always has an explicit temporal symmetry, and, finally, describes a genuine Lindblad evolution- 
even in the infinite-dimensional/continuous spectrum case. Application of the latter to a prototypical case of interest 
will be discussed, showing up, in contrast to the above mentioned pathologies of the CM approach, a reliable treatment 
of energy-dissipation and dephasing processes in semiconductor quantum devices. 

Contrary to standard master-equation formulations, ^-'^'^^ our approach leads to equations that are always of Lindblad 
typeji^ so that positivity is native in our adiabatic-decoupling strategy. Moreover, our approximation scheme holds 
true under the same validity regime of the conventional Markov approach: the so-called weak-coupling limit,— where 
the subsystem density matrix in the interaction frame moves slowly with respect to perturbative effects. 

The present article is organized as follows. After a general introduction, in Section U the key features of the formal 
CM approach and of the proposed alternative treatment are recalled and compared. Section|TT]presents the application 
of both schemes to a prototypical system, with the purpose of analyzing the failure of the CM and the success of the 
novel approach. The main conclusions, together with a brief summary, are contained in Section [IIII 



I. GENERAL FRAMEWORK 



To be more precise about the main features of the problem, it could be useful to recall its general formulation 
based on the fully operatorial approach proposed in Ref. [l^. Given a generic physical observable A — described 

by the operator A — its quantum plus statistical average value is given hy A ~ tv^Ap^, where p is the so-called 
density-matrix operator. Its time evolution is dictated by the global (system plus environment) Hamiltonian, that can 
be regarded as the sum of a noninteracting contribution, plus a system-environment coupling term: H — Ho + H'; the 
corresponding equation of motion for the density-matrix operator — also known as Liouville-von Neumann equation — 
in the interaction picture is given by: 

^it)^-^[n\t),p^ , (1) 

where denotes the interaction Hamiltonian H' written in units of h. 

The key idea, common to any perturbation approach, is that the effect of H' is "small" compared to the free 
evolution dictated by Ho- Following this spirit, by formally integrating Eq. ([1} from to to the current time t, and 
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inserting such formal solution for /3*(t) on the right-hand side of Eq. ([T]), we obtain an integro-differential equation of 
the form: 

^(^) = -^[7i^(^),/5X^o)]-^rf^'[?i'(^),[?i^(^'),/5X^')]] • (2) 

So far, no approximation has been introduced: Eq.s ([1]) and ([2]) are fully equivalent, we have just isolated the first-order 
contribution from the exact time evolution in Eq. ([T]). 

The final goal within the Markovian picture is an equation in which the time evolution of the complete density 
matrix, in Schrodinger representation and omitting the coherent rotation due to the first order contribution, can be 
expressed by means of a super-operator L^'"**"' as 

^ =L9'°'"^'(p). (3) 

In particular, by denoting with {|A)} the eigenstates of the noninteracting Hamiltonian Ho, and neglecting energy- 
renormalization contributions,— the effective equation ([3]) written in such a basis has the form: 

" 2 ^ ["^AiAa^A'iAiPA'iA^ - T'ai A^, ^A^ A; PA^ A2] + ^-C (4) 
A'jAj 

Equation ([4]) describes the evolution of the global -system {S) plus environment {E)- density matrix. However, 
in the study of electronic quantum phenomena in semiconductor nanostructures most of the physical quantities of 
interest depend on the electronic system only; the corresponding operators have therefore the form O — (8) 1^- 
Their expectation values may then be written as 

(O) = tr {Op} = tr I (6^ ® = tr {o^p^} (5) 

where 

- tr^ {p} , (6) 

is the reduced or electronic density matrix, obtained by tracing over the environment variables. 

The above mentioned reduction scheme can be interpreted in terms of a projection super-operator P, that acts on 
the generic global density matrix performing the partial trace over the environment degrees of freedom: 

Pp = tTE {p} (^p^ , (7) 

p^ being the environment density-matrix operator. This step is crucial to study the weak-coupling limit; as we shall 
see, it needs to be considered explicitly only later on, when talking about positivity. For the moment, we limit ourselves 

to assume that P (^H' {t),f'p\j ~ 0: This is physically justified for a wide family of perturbing Hamiltoniansjii and 

allows us to neglect the first order term on the right hand side of Eq. that gives no contribution when projectedi^ 
Hence for ease of exposition we shall now study all the projection-independent features of our model, safely neglecting 
first order terms. 

Since p^ is the only quantity entering the evaluation of the average value in Eq. ([S]), it is desirable to derive a 
corresponding equation of motion for the reduced density-matrix operator. However, P does not commute with the 
scattering super-operator one wants to project, thus hampering us in obtaining a closed equation of motion for p^ . 
Once at this point, to overcome such a problem, the typical procedure is to assume that the dynamics of the electronic 
subsystem does not significantly perturb the environment, which therefore is described by a time-independent density 
matrix. This can be achieved when the latter is, for example, maintained in a quasi-equilibrium regime and/or driven 
by some external mechanism. 

Let us now apply this approach to Eq. in the framework of Eq. ([7]): 

At dt 
Rewriting the latter on the basis given by the unperturbed electronic states {|q;)}, one obtains 

dn^ 

^Faia2 _ \ ^ u scatt S (q\ 
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where the complete super-operator 
is made up of four terms, that is 



^aia2,a[a'2 ~ 2\"'-"^'"'i"'2 aia\ ,0303 ^a'^a^ + 1 • (H) 

\ "3 / 



Here, the generahzed scattering rates V are obtained applying the projection P to Eq. (|4]). In this respect, diverse 
approaches may be employed. In particular, in the following Subsection II Bl the generalized scattering rates V will 
derived according to an alternative procedure for the Markov limit, able to overcome the well-known critical issues of 
the conventional one. To better compare the key features of both schemes. Subsection II Al contains a brief description 
of the latter. 



A. Conventional Markov limit 



Let us therefore focus on the time integral in Eq. ([2]). Here, the two quantities to be integrated over t' are 
the interaction Hamiltonian T-C and the density-matrix operator p*. In the spirit of the perturbation approach 
previously recalled, the time variation of can be considered adiabatically slow compared to that of the Hamiltonian 
i-L written in the interaction picture, i.e., W{t') = tjl{t')'Hij^{t'); indeed, the latter exhibits rapid oscillations due 

to the noninteracting evolution operator Uo{t) — e Therefore, in the standard (and problematic) Markov 

approximation the density-matrix operator /5' in the interaction frame is simply taken out of the time integral and 
evaluated at the current time t. 

Following such a prescription, the system dynamics written in the Schrodinger picture for the case of a time- 
independent interaction Hamiltonian % comes out to be: 



dt =-^[^'^1-2 



with 



(12) 



(13) 



t-to 



The effective equation in (|12l) has still the double-commutator structure in ([3]) but it is now local in time. In this 
scheme, the generalized scattering rates of the Eq. ([3]), obtained within the completed-coUision limit to — > — oo, turn 
out to be 



27r 



(14) 



where e\ denotes the energy corresponding to state |A). 

Once projected with P, and considering a density matrix in the form of Eq. ([7]), Eq. 



T2]) reduces to 
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z[Ho,/5^]--P([H',[/C,p^]]). 



(15) 



The Markov limit recalled so far leads to significant modifications in the system dynamics: while the exact quantum- 
mechanical evolution in Eq. ([T|) corresponds to a fully reversible and isoentropic unitary transformation, the instan- 
taneous double-commutator structure in Eq. (jl2p describes, in general, a non- reversible (i.e., non unitary) dynamics 
characterized by energy dissipation and dephasing. However, since any effective Liouville super-operator should 
correctly describe the time evolution of p and since the latter, by definition, needs to be trace-invariant and positive- 
definite at any time, it is imperative to determine if the Markov super-operator in Eq. (|12p fulfills this two basic 
requirements. 

As far as the first issue is concerned, in view of its commutator structure, it is easy to show that this effective 
super-operator is indeed trace-preserving. In contrast, as discussed extensively in Ref. [19[, the latter does not ensure 
that for any initial condition the density-matrix operator will be positive-definite at any time. This is by far the most 
severe limitation of the CM approximation. 
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The well-known semiclassical or Boltzmann theorys can be easily derived from the quantum-transport formula- 
tion presented so far, by introducing the so-called diagonal or semiclassical approximation. The latter consists in 
neglecting all non-diagonal density-matrix elements (and therefore any quantum-mechanical phase coherence between 
the generic states Ai and A2), i.e., PA1A2 = /ai^AiA2j where the diagonal elements fx correspond to the semiclassical 
distribution function over our noninteracting basis states. Within such approximation scheme, the quantum-transport 
equation (|12p reduces to the well-known Boltzmann equation: 

^=E(^aa'/a'-Pa'a/a) , (16) 

A' 

where 

27]- 

Paa' = Paa.A'A' = — \H'^y\^5iex-ey) (17) 

are the conventional semiclassical scattering rates given by the well-known Fermi's golden ruleJ^ 

At this point it is crucial to stress that, contrary to the non-diagonal density- matrix description previously intro- 
duced, the Markov limit combined with the semiclassical or diagonal approximation ensures that at any time t our 
semiclassical distribution function f\ is always positive-definite. This explains the "robustness" of the Boltzmann 
transport equation in (|16p . and its extensive application in solid-state-device modeling as well as in many other fields, 
where quantum effects play a very minor role. In contrast, in order to investigate genuine quantum-mechanical phe- 
nomena, the corresponding CM super-operator L^'^J"' obtained from Eq. (jl2p cannot be employed, since it does not 
preserve the positive-definite character of the density matrix p\^\^. 

As highlighted before, the origin of these pathologies is in the lack of temporal symmetry in the construction of the 
CM approximation; namely, in the double commutator expression in Eq.s ()12|) and (|15|) . Motivated by this evidence, 
in the next Subsection, we shall propose an alternative adiabatic approximation method in which the hamiltonian 
operators of Eq. ([2]) are always time-balanced, so that the temporal symmetry remains always explicit and is never 
broken. 



B. Alternative approach to the Markov Hmit 



To introduce our alternative formulation of the problem, let us go back to the integro-difFcrcntial equation in 
We can formally integrate the latter from to to t, and obtain 



p\t)^p\U)-l / dil W{h),p\to) 



dti / dt2 

to 



(18) 



This equation, named after Sadao Nakajima and Robert Zwanzig, is still exact, and it is nothing but the Dyson 
equation after two iterations. 

As far as the first integral on the right-hand side of Eq. is concerned, the most it can cause is an energy 
renormalization effect that, after the partial-trace projection in our proposed scheme, will straightforwardly result 
into a coherent term. Since the latter does not affect positivity in any case, we shall forget about it in the following. 

Let us now perform a change of variables. In particular, starting from the two times ti and ^2, we introduce a 
"relative" time 



and a "macroscopic" time 



T 



to 



t2 



(19a) 



(19b) 



This change of variables has very sound bases, as it is typical and well established for a wide variety of contests, such 
as Wigne r's p hase-s p ace formulation of quantum mechanics?^ standard quantum kinetics Green's function methods 
(see e.g., [22 and|23|) and even classical radiation theory (e.g., in the treatment of Bremsstrahlung) . The basic idea 
is that the relevant time characterizing/describing our effective system evolution is the macroscopic time T. 

Having this profound difference between T and r in mind, we can rewrite Eq. (|18p as: 



9(t,T) 



JdT Jdr[w{T+l),[w{T-l),p'{T-l 



(20) 
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where the function g comes from the domain of integration in the (T, r) plane and can be expressed as^l 

9{t.T) = 



t 








T - 


(0 


2 ^ 







■ (21) 

The integral over the relative time r in Eq. (I^Ul) requires some manipulation. In particular, by introducing a proper 
gaussian cut-off factor, e 2t^, the upper limit g{t,T) may be extended to infinity; that is 



9{t.T) 



/(r)dr 



fir)e- 



dr. 



(22) 



In other words, the last integral on r gets effectively cut off at the collision time t, having in mind that the contributions 
resulting from the double commutator structure get negligible at longer times. In particular, t scales up as some 
negative power of the coupling constant (that is the ratio between the interaction matrix element and the unperturbed 
energy gap) and goes to infinity when the interaction is weak. 

In this framework, we now evaluate the first derivative of Eq. ()20|) with respect to t and obtain: 



dt 



dr e 



(23) 



In the spirit of the adiabatic approximation previously recalled, the density-matrix operator , being a slowly varying 
function of r, can be taken out of the time integral and evaluated at the current time t. 

The skew-adjoint part of Eq. (I23p -the so-called scattering part- is our main interest, the self-adjoint part being 
just an energy-renormalization term that does not threaten positivity; we then get: 



dr. 



(24) 



It is evident how the proposed time symmetrization gives rise to a fully symmetric super-operator, contrarily to the 
strongly asymmetric Markov super-operator in Ref. 18"]. 

The second crucial step to get a genuine Lindblad super-operator for the global dynamics is to exploit once again 
the slowly- varying character of the density- matrix operator on the right-hand side of Eq. (j24p . The key idea 
is to perform on both sides of Eq. (f24| a so-called temporal "coarse graining", i.e., a weighted time average on a 
microscopic scale where the variations of (t) are negligible. In the weak and intermediate coupling regime this time- 
scale is comparable to t; we shall then perform the above mentioned convolution employing once again a gaussian 
correlation function of width |, i.e.. 



dr' 



'27rt 



dr e 



W it 



(25) 



Moving back to the original Schrodinger picture and combining the two gaussian distributions, the above equation 
can be rewritten^^ in the following compact form: 



witb2^ 



dp 
dt 



(26) 



C = 



/27rt 



dt' Wit') 



(27) 



This is the genuine Lindblad-like super-operator we were looking for; indeed, the operators £ are always Hermitian, 
and therefore the effective dynamics they induce is positive-definite. 

Let us finally rewrite the new Markov super-operator in Eq. (|26p in our noninteracting basis {|A)}, defined by the 
(possibly generalized) eigenvectors of Ho . In this way, we now obtain an effective equation of motion in the same form 
of Eq. Q with symmetrized quantum scattering rates 



—H' H'* 

^ AiAj A2A2 



27re 



exp ■ 



(cAi 



(cas - ^a; 



4e^ 



(28) 
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Here, e = fi/t = fiZU is a measure of the energy uncertainty in the interaction process induced by our temporal coarse 
graining. 

The above scattering super-operator can be regarded as a generahzation of the conventional Fermi's golden rule to 
the density matrix formalism; Indeed, in the semiclassical diagonal case (Ai — A2, A'j^ = Aj) it reduces to what could 
be considered a dressed vertex-smoothed version of the Fermi's golden rule 

Pav =Paa,va' = -i^e-(— v)V2.^ (29) 

n v27re 

In the limit of infinite correlation-time (w — > 0), the standard scattering rates given by Eq. (jl7p are readily recovered. 

At this point, one may wonder whether the approximate dynamics so obtained suffers from lack of generality, the 
gaussian smoothing and cut-off having been put by hand. However, it is important to stress that one has to limit 
the analysis to the asymptotic features, as the approximation is valid in the weak-coupling regime. As such, and 
without loss of generality, our gaussian choice, among all the possible asymptotic markovian approximations of the 
exact dynamics that guarantee a positive evolution, is a good representative. Moreover the use of gaussian functions 
has the marked advantage of allowing for the analytical treatment we are carrying out in the present article. 

The transition rates in Eq. (j28p could be regarded, in some sense, as a quantum version of the celebrated Fermi's 
golden rule. This should not generate confusion: of course the scattering rates obtained with the latter are intrinsically 
quantum in their derivation, but, once computed, they give rise to the Boltzmann equation in (|16l) . which describes a 
classical Markov process^^ for classical probabilities; the evolution of the probability density function f\. Conversely, 
the transition rates in Eq. (1281) do not describe a classical Markov process, but rather its quantum analog: a so called 
quantum dynamical semigroup^ for the full density matrix. In other words, the dynamics of the entire (i. e., diagonal 
plus non-diagonal terms) reduced density matrix is determined by a Lindblad-like closed set of equations of the form 

^ = -z[Ho,/5^]-ip([A[Ap^]]), (30) 
and it is still completely positive and trace invariant. 



II. APPLICATION TO A PROTOTYPICAL SYSTEM 



To highlight the critical issues in the CM approximation, let us consider as system S a semiconductor quantum dot 
(QD) interacting with a solid-state environment E and forced by an external single sound wave source of wavevector 
Qo . The free Hamiltonian for such a system has the form 

Ho=Hs®'^E + is®HE = Y^ hojacicc, KK' (31) 

a q 

In Eq. ([3T|). the bosonic operator b^^ (5q) denotes creation (destruction) of a phonon excitation with wavevector q and 
energy Cq = hujq, while the fermionic operator cj^ (cq,) denotes creation (destruction) of an electron in the state a 
with energy = fiuJa- 

The noninteracting (carrier-plus-quasiparticle) basis states |A) are given by the tensor product of electronic states 
I a) and quasiparticle states K'T.q}) corresponding to the occupation numbers {uq} 

|A) = \a) ® |{nq}) ; (32) 

the noninteracting energy spectrum is then the sum of the electronic and quasiparticle energies 

hldX = huJa + f^q nq . (33) 

q 

The perturbation Hamiltonian H' can be written as: 

H' = hJ2 {Hqbq + Hljfelj) = /i(H"'' + ?i"") , (34) 
q 

here, Tiq ~ Ti^-q are electronic operators (parametrized by the quasiparticle wavevector q) acting on the a subsys- 
tem only. The two terms in Eq. p4p -corresponding to quasiparticle destruction and creation- describe electronic 
absorption {'H"''') and emission {H'^"^) processes. 
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In this discussion we write 'H' as: 

'^'=J2Y1 5aia2,q ® 6q + ll.C. = (sai ,q + Sa^ai ,q ^q) , (35) 

Qia2 q OLia2 q 

that is, Tiq = J2aia2 9°'io'2,<i cJ,jCq,2, for a generic complex-valued coupling coefficient gaia2.cii depending on the 
wavevector q. In particular, within the framework of a deformation potential interaction, the following form for the 
5aia2,q applies: 



W2,q = ^27^4j; / C,(r)e^''-''^a.(r)dr = 5q(ai|e'^i-|a2) , (36) 

where e is the deformation potential, p is the crystal mass density, and 0Qj(r) [(/'a2('')] the wavefunction of the QD 
single electron state ai [a2]- 

The coherent contribution can be taken into account by formally introducing a free-part rotation, deriving from 
Ho, so that L = L''<='= + L"'^^". To write/have the super-operator in matrix form, we introduce the matrix-to-vector 
mapping i— > g^' given by 

£'f — ^ Pfia2 '^i.(Ql-l)Tl + Q2 ■ (37) 

aia2 

The dynamical equation for the reduced density matrix can then be rewritten as 

^^g^''{t)^hg^{t). (38) 

In the estimation of L^*^***', we keep in mind that the projection on the environment degrees of freedom of the 
operators fojj6q' and &q'6q has the following, finite, value 

trB{6^6q'/5^} =iVq(5qq,, (39a) 

tiE {K-Kp^} = (^q + l) -^qq' ' (39b) 

with A^q denoting the average occupation number for state q. Moreover, coherent phonon contributions are not 
accounted for, that is 

tr^jSq/)^} =tr£{6tpS| ^0 (4Q) 

For the sake of simplicity, let us now assume the effect of the external sound wave source to be dominant in the 
wavevector distribution of the phonon population. This is indeed a reasonable assumption when, e. g., operating at 
low lattice temperature. Such a regime may be effectively described by a time-independent diagonal density matrix 
for the environment, having one element overwhelming the others. It is important to stress that this does not mean a 
strict one-phonon environment (where any Markov approach would be meaningless), but, on the contrary, a situation 
in which, while the electron system is forced to very efficiently couple to a certain phonon mode, the latter interacts 
with the whole lattice. 

In this respect, from now on, we limit our attention to the interaction between the system and the externally induced 
sound wave only. That is, the summations over q will be evaluated just for the phonon-modc qo. In particular, this 
amounts to say that A'q ~ (5(q — qo). 

A. Failure of the conventional Markov approach 

In the CM limit, the super-operator L in Eq. © turns out to be 

L"'^^" , , - -27r V f^^^^ 6 S A+2n(A^^'' +.4^^^ (Al) 

^aiQ2,aia^ ^" a'^a^, 01302 ^la^ ^ -^0103, aaa'i "202^ ' ^" \^"^Q^Q2,aiQi ^ -^0^02,010'^^ ' V^^^ 

as 

where 
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and 



(42b) 



In the latter equations / are defined as 



(q) = 3a2«i(q)5a;a;,(q), 



(43) 



and L0aia2 — '-^ai ~ ■ 

Indeed, using the properties of the Dirac-(5 function 



- go) 

|VL.(qo)| ' 



(44) 



where lo{(\o) = uJo — T'^aiQ2i can easily perform the integral in Eq. ()42p . Moreover, we may consider our acoustic 
wave qo in the long-wavelength limit so that we can approximate 



w(q) w |Vw(q)| (?, 
to have an explicit expression for the typical integrals A: 



A 



(1) ^ ^|q° 



E 



Ql Q2 iCti CKq 



1 1 

2 + 2 



(45) 



(46a) 



and 



(2) _ 



E 



' ai a2 ,CKi cto 



1 1 

2 + 2 



(46b) 



where the functions / (or similarly g) have to be evaluated in qo, and the Kronecker— (5 accounts for energy conser- 
vation. 

Now, to further focus on our prototypical case of interest, we limit our discussion to the bound-to-continuum 
transitions induced by the external source in the electron system of a QD having a single confined state. In this way, 
labelling the latter as |1), the only allowed transitions are those into a final state |2) so that w(qo) = lOo = ^21- 

The complete super-operator L = L^'''"' + L'^'^'^** may then be written as: 



L = 





iw2i 

-iuj2i 

GO 



47ra;o 



c 



2^ 





2N+1 
2 

„* 2jV+l 
'/ 2 







,2^ 

_ 2JV+1 
2 





(7V + 1) 

^ _2 
f-* jV+l 

~(iV + l)/ 



(47) 



where r] 



gi2(qo) 
321 (qo) 



gii(qo)-g22(qo) 
921 (qo) 



, and the dimensionless ratio cr, defined as 



V 

8^ 



l32i(qo)|' 



(48) 



measures the relative strength of the perturbation Hamiltonian. 

As can be seen by explicit inspection of the matrix elements of L, the equations of motion for the populations 
and P22 ^i'^ completely decoupled from those referring to the polarization elements p^i and pfi — {Pi^T ■ the 
contrary, the presence of the C and t] terms lead to a polarization evolution that strongly depends on the populations 
and can be divergent in some cases. 

Defining, for later convenience, the quantities 



K = 27r(27V+ l)cr 

At' = |?7|^K, 



(49a) 
(49b) 
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we can write the four eigenvalues of L as: 

Ml = 

^2 = -2a;o K (50) 

^3^4 — UJo ( — K ± Vkk' — l) . 

The null eigenvalue, /ii, guarantees the occurrence of a steady-state solution. It can be shown that the latter corre- 
sponds to the vector 

Q^^^ (iV+l,0,0,iV)^, (51) 

2iV-f 1 

that is, to the 2x2 matrix 

r - I I . (52) 

We will later show that, in spite of the fact that such a steady state seems so physically sound, the latter may not 
ever be reached. As a consequence, the already known pathologies of the CM approach, such as the possible lack of 
positivity,^ will appear to be much more serious than commonly believed, as they will in fact be proved to show up 
not only in some initial transient, but also, and heavily, in the long-time/steady state regime. 

Coming now to the remaining eigenvalues, by inspection of Eq. (1501) it is evident that the behaviour of /ia (corre- 
sponding to the -t- sign) is crucial. In particular, when 




Vkk' - 1 > k , (53) 

5R(/i3) gets positive and this may lead to a divergence in the polarization terms. After straightforward manipulation, 
Eq. ((55)) may be cast into the following form 

k2(|,^|2 _ 1) > 1. (54) 

The latter is surely not verified if \rj\ < 1, a case in which the steady state in Eq. ([52]) is always eventually reached. 
However, when > 1, by rewriting Eq. (|54l) with use of Eq.s (|49l) . one gets 

11- 

a>— = a . (55) 

2N + 1 27r^|77|2 - 1 

that is, 3?(/i3) gets positive for a larger than a critical value a. 

An alternative approach to Eq. (|55p -and more consistent with the fact that external source intensity is an experi- 
mentally controlled degree of freedom- is to read it in terms of N: 

iV>i(i^^^_-l) ^iV. (56) 



2 U 2ti^\ 



Figure [T] shows the real part of the four eigenvalues /i^ of the super-operator L given by Eq. (|50|). in units of Wo, 
as a function of <t/ct, for the case rf' = 1.31 (this choice will be motivated later on), that is \rj\ > 1. Besides the 
null eigenvalue, fXi, and the negative-slope one, /X2, the detrimental behaviour of 9? (/X3), discussed above, is evident. 
To have a deeper insight into the critical aspects that this may cause, let us now further specify the features of our 
prototypical system. 

To provide a quantitative analysis, the QD three-dimensional confinement potential V is assumed to be made up of 
a finite square-well profile V", of width d, along the growth (z) direction, and a truncated two-dimensional parabolic 
potential in the xy plane. This is a reasonable configuration for state-of-the-art Ill-V-based self-aggregated QDs. 

For the sake of simplicity, and without compromising the analysis we are interested in, in the following we will assume 
that the factorized form, typical of infinite confinement, remains valid also for our truncated QD. In particular, as 
single particle bound state |1), we shall consider the lowest energy one (r|l) = V'oi(r) — 4>^{x, y) (2), where (j)\{z) is 
the ground state of the one-dimensional quantum- well along the z-axis, centered in z = and of width d, and (j)^ (x, y) 
is the ground state of the harmonic oscillator in the xy plane. Finally, the localization energy (measured from the 
barrier conduction band minimum) E\oc — ^t^ioc of state V'oi(r) is set to 7 meV. 
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Figure 1. Real part of the eigenvalues of the super-operator L, in units of cjo, as a function of the dimensionless parameter ^ 
for the CM approximation. Data refer to the case |77p = 1.31 and TV = 9, leading to a = 0.01. 



To compute the entries of the matrix L, namely 77, C and cr, we have to calculate the form factors gaa' (q) coupling 



ipoi{r) to the continuum state (r|2) = e /vk. In particular, fixing the wave-vectors qo and k along the growth 
direction (z), the quantities 77, ^ and a turn out to be 

coa[d{k + qo)/2] £■ {k - qof - tt'^ 



V 



C 



cos[d{k 



qo)/2] 



(P{k + qoY-Tr^ 

qo)"^) qo sec [d{k — qo)/2\ sin [dqo/2\ 



27r3/2 (47r2 - d^ql) 
cos[d{k — (?o)/2] 



7r2 - d^{k - qof 



(57a) 
(57b) 
(57c) 



where a — [m* u) / 2h) ^ 



being the paraboloid curvature and m* the electron effective mass. 



Let us now consider an initial condition in which the QD bound level is occupied. This configuration, in the {|1), |2)} 
basis, is represented by the following density matrix 



1 




(58) 



In the assumed low-temperature limit, and in the absence of any external pumping source, a QD prepared in such 
a state would indefinitely remain in the latter. This is not the case when the sound wave source is on. At this 
point, however, some severe pathologies of the conventional approch may show up. Indeed, once the density matrix 
in Eq. (j58|) is rewritten in the basis given by the eigenvectors of L, we do expect that the contribution arising from 
the component p'^-^'^^) may lead to a divergent dynamics when the critical conditions of Eq. (j53p are met. 

A convenient way to picture out such a non trivial behaviour is given by Fig. [2j In the latter, N is plotted as a 
function of the phonon energy huJo- This produces a partition of the {N , fnoo) parameter plane into distinct regions, 
according to whether the condition in Eq. (j56|) is satisfied or not. In particular, there appear a "nonphysical region" 
(marked in grey), where the density matrix eigenvalues diverge, and a region in which, on the contrary, Eq. (1561) is 
not satisfied therefore guaranteeing the occurrence of the expected steady-state solution. 

To better focus on such a density matrix dynamics, let us consider two diverse setups, differing in the energy, Swo, 
of the incident phonon beam: a first one in which hwo — 11 meV, and a second one with fiwo= 10 meV. For the 
former, marked by a diamond in Fig. [2J the critical region shows up for N > 10; this is not the case of the latter, 
marked by a triangle, which falls in the 'safe' region. Indeed, despite of the fact that the two phonon energies are 
quite similar, the effects on the system evolution are dramatically different, as we are showing in the following. 

In particular, considering the time evolution of the density matrix elements p^- when our electronic system is initially 
prepared in the Pg state and /ia;o= 10 meV, the population terms, pfi and P22J ^'I'S found to converge to a steady-state 
value in a sub-picosecond timescale, while the polarization contributions, and pfj^, decay exponentially to zero, as 
expected. 

Let us now move to the case in which fiWo = 11 meV. As we may expect from the evidences given in the previous 
Section, some critical aspects should appear in the dynamics of p^ making the latter no longer positive definite. Indeed, 
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Figure 2. Partition of the (A^, huio) parameter space produced by the N = N equation. The triangle (diamond) highhghts the 
configuration in which the phonon energy fkuo amounts to 10 meV (11 meV), with A'' = 20. The circle corresponds to the case 
&Jo = 11 meV and iV = 9 (a- = 0.01). 

as shown in Fig.[3l though the population time evolution (continuous lines) still appears somehow reasonable (i.e., the 
trace is preserved), the polarization one (dashed and dot-dashed lines) is unphysically divergent. As highlighted before. 
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Figure 3. Time evolution of the system density matrix elements, pfj, in the {|1), |2)} basis, when the QD is initially prepared 
in the |1) state. Continuous lines refer to the populations pfi and P22', the dashed (dot-dashed) line corresponds to the real 
(imaginary) part of the polarization term pi2. The phonon frequency is huio = 11 meV, a = 0.0092, and TV = 20. 

the origin of these pathologies is in the lack of temporal symmetry in the construction of the CM approximation; 
namely, in the double commutator expression in Eq. (|12p . 



B. Success of the alternative Markov procedure 



We are now interested in comparing the results of Subsection lll Al with those derived from above described alternative 
approach to the Markov approximation; we will show that the previously mentioned pathologies of CM are now 
brilliantly overcome. 

First of all, performing the trace over the environment degrees of freedom, Eq. (|30p turns into 



(59) 
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with 



where 



(60) 



(61) 



5Qia2,q = ffaiQ2,q> 5<tiQ2,q = 3a2Qi,q' ^^d is again the environment density matrix. 

With the above exphcit expression of the interaction Hamihonian, it is easy to perform the time integrals in Eq. (|60p , 
since the latter only involve products of gaussian functions and plane waves. Equation (j59p is then cast into the form 
of Eq. ([9|), however, contrary to the CM result, the novel V terms are now given by^S 



°2°2/ 



q:± 



(62) 



In the completed-collision limit (lj 0), that is, when the broadening due to the the gaussian correlation function 
is null, the first exponential term in Eq. (|62p results into a Kronecker (5-function, while the second one produces a 
Dirac ^-function as follows 



^2°2 87r3 



^ /dq TVq + - ± - (q) (q) 5 



(63) 



While the former forces the system to have a discrete spectrum -the only one considered in the exiting literature, 
up to now- the latter corresponds to energy conservation. It is straightforward to verify that the usual semiclassical 
rates Vaa' given by the Fermi's golden rule may be obtained from Eq. (155)) when ai — a2 — X and a'^ — a'2 — \' ■ 

For a more quantitative analysis, let us now focus on the prototypical system considered in Section |lTl that is a 
semiconductor QD interacting with a collimated and monoenergetic source of phonons of wavevector qo. The latter, 
in particular, only couples the QD bound state |1) to a specific continuum one, |2). Substituting for the explicit 
expression TVq = iVg^ J(q — qo), the generalized scattering rates become 



(64) 



where a = 1,2. In our approach, the markovian super-operator L — Lfioc _j_ j^scatt j^g^y .^^^ written in matrix form. 
In particular, the free part, L^'''''^, is identical to the one obtained in Subsection III A[ 



free 



/o 00 

iuj2i 
-iuj2i 

Vo 



(65) 



while L*" 



now becomes 



( -N 



V N 





2N+1 
2 










2N+1 
1 

-N 



A^ + 1 \ 




(66) 



with a as given in Eq. pS)) . 

Let us now compare Eq. ([M]) with the scattering part of Eq. (|T7| in the CM-limit case. First of all, when |5r2i| — 0, 
the scattering super-operator vanishes in both of them. However, in Eq. (j66|) this result simply derives from the 
fact that L^'^'i'* depends only on (721: through a. The lack of any coupling between the populations and polarizations 
dynamics in Eq. (j66p leads to a time evolution of p that, contrarily to the CM approximation, does not have any 
contribution^^ Moreover, although the first and last rows of the L'^'^'^** matrix are identical in the two approaches. 
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the dynamics of the polarizations is totally different. This feature is clearly evident when one diagonalizes L. Indeed, 
the eigenvalues of the latter are 





= 








= imT 


^4 





(67) 



where k is related to tr by k = 2tt{2N + l)a. 

Figure E] presents the eigenvalues /ii just obtained as a function of a /a. This choice is made to directly compare 
these results with the ones showed in Fig. [T] in particular, a has been defined in Eq. (|55p . within the CM-limit 
framework, and is here set to the value a — 0.01. We stress that the real part of the eigenvalues in Eq. (l67l) is always 
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Figure 4. Real parts of the eigenvalues fii of the super-operator L, plotted as a function oi a/a to allow for a direct comparison 
with Fig. [T] 



non-positive, regardless of the value of the gij's; the zero eigenvalue corresponds to the steady state solution 



1 



2N + 1 



N + 1 



N 



(68) 



identical to the result in Eq. ([5^ . obtained with the CM approach. 

Let us now turn to Eq. (l38l) . By inspection of Eq.s (|65p and (l66l) L contains an inner 2x2 diagonal block, relative 
to the polarizations dynamics, while the evolution of the population is dictated by coefficients identical to those of 
the CM approach. This strongly simplifies the problem, allowing for an analytical solution. The vector is then 



9' it) 



^•S,(mi 

\ P22 



(pfi(O) 
(pfi(O) 



~1k uJot 



Pii 

Pi2(0) e-(''+*)'^°* 



(69) 



where pf, (0) are the elements of the 2x2 density matrix at time t = 0. It is important to notice that can now 

and a contribution Sp{t) which 

2ri, 

(70) 



be read as a sum of the steady-state, i.e., t-independent, solution given in Eq. 
decays exponentially and contains the characteristic times Ti = (2k Wo)"^ and T2 

p^t) ^p' + Sp{t) . 

In other words, whatever being the initial condition p^{0), the density matrix reaches the equilibrium when the 
real parts of the exponential functions in Eq. (|69p approach zero; the intensity of the 6p{t) term is proportional to 
the difference between the initial condition and the steady state solution. Moreover, contrarily to the CM case, where 
severe instabilities may result from a diagonal perturbation produced in the system steady-state, here the form of 
6p(t) is so that the latter does not affect the polarization dynamics at any time. 

To be more quantitative, and also to compare the results of Subsection III Al with the dynamics given by our 
alternative approach, let us consider a density matrix p^ initially prepared as in Eq. (j58p and the configuration 
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marked by a diamond in Fig. The latter was a critical situation for the CM method that produced the divergent 
behaviour in the polarization dynamics shown in Fig. [31 The situation resulting from our approach is completely 
different, being stable and physically meaningful. In particular, the time evolution of given by Eq.s (|65|l and (|66|) 
produces the results shown in Fig. [5] The population terms pf^ and pfj decay exponentially from the initial values 1 
and 0, respectively, to the steady state ones pfj^ ~ 0.51 and pf2 — 0.49, which are identical to those found in the case 
reported in Fig. [3] for the CM approach. The main discrepancy between the two methods appears when considering 
the behaviour of the polarization terms: following our prescriptions, and conversely to the divergent CM result, the 
latter remain null, as they initially are. 
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0.2 
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0.00 0.05 0.10 0.15 0.20 

t(ps) 

Figure 5. Time evolution of the population and polarization terms of , in the |2)} basis, when the QD is initially 
prepared in the |1) state. Continuous lines refer to the populations pfi and p22', the dashed (dot-dashed) line corresponds to 
the real (imaginary) part of the polarization term pf2- The phonon frequency is huo = 11 meV, a — 0.0092, and N = 20. 



III. SUMMARY AND CONCLUSIONS 

We have presented and discussed a general density-matrix description of energy-dissipation and decoherence phe- 
nomena in open quantum systems, able to overcome the intrinsic limitations of the conventional Markov approxi- 
mation. More specifically, contrary to the usual single-particle correlation expansion, we have investigated the effect 
of the Markov limit, in a fully operatorial approach, before applying any reduction procedure. This has allowed us 
to better identify the general properties of the scattering super-operators entering our effective quantum-transport 
theory and to propose an alternative adiabatic scheme that does not threaten positivity at any time. The key idea 
of our approach consists in performing the temporal symmetrization and coarse graining of the scattering term in 
the Liouville-von Neumann equation, prior to the reduction over the environment degrees of freedom. The result is 
a robust treatment of energy-dissipation and dephasing phenomena in state-of-the-art semiconductor quantum de- 
vices, able to describe a genuine Lindblad-like evolution that also recovers the Fermi's golden rule features in the 
semiclassical limit. 

Applications to the prototypical case of a semiconductor quantum dot exposed to a single-phonon source are 
discussed, highlighting the success of our formalism with respect to the critical issues of the conventional Markov 
limit. In particular, we have shown very important intrinsic limitations of the CM approach when employed for large 
time and/or steady state analysis, producing severe instabilities that show up in divergent polarization terms. 

In the present article we focused on the dissipative part of the scattering super-operator, deliberately neglecting the 
energy level renormalization. Indeed, as can be checked both analytically and numerically, the effect of the latter is 
to worsen the pathologies of the conventional Markov limit, rather than to remedy for them. Moreover, the presence 
of energy renormalization effects in our proposed scheme would result into a coherent term in the dynamical equation 
for the global density matrix, which does not affect positivity in any case. 

We stress that our formulation significantly generalizes preexisting theories, as it gives a positive dynamics for a 
considerably large class of projections (i.e., ways to chose the subsystem), irrespective of the subsystem's dimension or 
spectral properties. This allows, on one side, to investigate subsystems with both discrete and continuous spectra -a 
feature largely shared by mesoscale electronic and opto-electronic quantum devices- and suggests, on the other side, 
a new way to treat electrical contacts for quantum devices, thus revealing an interesting potential for an alternative 
quantum transport formalism. 
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Finally, it is imperative to stress that, in the presence of a strong system-environment interaction, the adiabatic 
decoupling investigated so far needs to be replaced by more realistic treatments, expressed via non-Markovian integro- 
differential equations of motion^ (i.e., with "memory effects"). Again, while for purely atomic and/or photonic systems 
it is possible to identify effective non-Markovian evolution operators j'^^ for solid-state quantum devices this is still an 
open problem. 
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